function [bestLog2c,bestLog2g,bestcv] = autogrid(trainLabel,trainData,options)
% auto find the best C and gamma for RBF kernel
% can define more param options in `options`

if(nargin < 3)
    options = '';
end

findMax = true;
bestcv = 0;
if (strfind(options,'-s '))
    findMax = false;
    bestcv = 999999;
end

stepSize = 1;
log2c_list = -20:stepSize:15;
log2g_list = -20:stepSize:20;

numLog2c = length(log2c_list);
numLog2g = length(log2g_list);
cvMatrix = zeros(numLog2c,numLog2g);

for i = 1:numLog2c
    log2c = log2c_list(i);
    for j = 1:numLog2g
        log2g = log2g_list(j);
        % -v 5 --> 5-fold cross validation
        param = ['-q -v 5 -c ', num2str(2^log2c), ' -g ', num2str(2^log2g), options ];
        cv = svmtrain2(trainLabel, trainData, param);
        cvMatrix(i,j) = cv;
        if(findMax)
            if (cv >= bestcv),
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end                
        else
            if (cv <= bestcv),
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end
        end
        
        % fprintf('%g %g %g (best c=%g, g=%g, rate=%g)\n', log2c, log2g, cv, bestc, bestg, bestcv);
    end
end

disp(['CV scale1: best log2c:',num2str(bestLog2c),' best log2g:',num2str(bestLog2g),' accuracy/rmse:',num2str(bestcv),'%']);

% Plot the results
figure;
imagesc(cvMatrix); colormap('jet'); colorbar;
set(gca,'XTick',1:numLog2g)
set(gca,'XTickLabel',sprintf('%3.1f|',log2g_list))
xlabel('Log_2\gamma');
set(gca,'YTick',1:numLog2c)
set(gca,'YTickLabel',sprintf('%3.1f|',log2c_list))
ylabel('Log_2c');


% ###################################################################
% cross validation scale 2
% This is the medium scale
% ###################################################################
prevStepSize = stepSize;
stepSize = prevStepSize/2;
log2c_list = bestLog2c-prevStepSize:stepSize:bestLog2c+prevStepSize;
log2g_list = bestLog2g-prevStepSize:stepSize:bestLog2g+prevStepSize;

numLog2c = length(log2c_list);
numLog2g = length(log2g_list);
cvMatrix = zeros(numLog2c,numLog2g);
bestcv = 0;
for i = 1:numLog2c
    log2c = log2c_list(i);
    for j = 1:numLog2g
        log2g = log2g_list(j);
        % -v 3 --> 3-fold cross validation
        param = ['-q -v 5 -c ', num2str(2^log2c), ' -g ', num2str(2^log2g), options];
        cv = svmtrain2(trainLabel, trainData, param);
        cvMatrix(i,j) = cv;
        if(findMax)
            if (cv >= bestcv)
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end                
        else
            if (cv <= bestcv)
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end
        end
        % fprintf('%g %g %g (best c=%g, g=%g, rate=%g)\n', log2c, log2g, cv, bestc, bestg, bestcv);
    end
end

disp(['CV scale2: best log2c:',num2str(bestLog2c),' best log2g:',num2str(bestLog2g),' accuracy/rmse:',num2str(bestcv),'%']);

% Plot the results
figure;
imagesc(cvMatrix); colormap('jet'); colorbar;
set(gca,'XTick',1:numLog2g)
set(gca,'XTickLabel',sprintf('%3.1f|',log2g_list))
xlabel('Log_2\gamma');
set(gca,'YTick',1:numLog2c)
set(gca,'YTickLabel',sprintf('%3.1f|',log2c_list))
ylabel('Log_2c');



% ###################################################################
% cross validation scale 3
% This is the small scale
% ###################################################################
prevStepSize = stepSize;
stepSize = prevStepSize/2;
log2c_list = bestLog2c-prevStepSize:stepSize:bestLog2c+prevStepSize;
log2g_list = bestLog2g-prevStepSize:stepSize:bestLog2g+prevStepSize;

numLog2c = length(log2c_list);
numLog2g = length(log2g_list);
cvMatrix = zeros(numLog2c,numLog2g);
bestcv = 0;
for i = 1:numLog2c
    log2c = log2c_list(i);
    for j = 1:numLog2g
        log2g = log2g_list(j);
        % -v 3 --> 3-fold cross validation
        param = ['-q -v 5 -c ', num2str(2^log2c), ' -g ', num2str(2^log2g), options];
        cv = svmtrain2(trainLabel, trainData, param);
        cvMatrix(i,j) = cv;
        if(findMax)
            if (cv >= bestcv)
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end                
        else
            if (cv <= bestcv)
                bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
            end
        end
        % fprintf('%g %g %g (best c=%g, g=%g, rate=%g)\n', log2c, log2g, cv, bestc, bestg, bestcv);
    end
end

disp(['CV scale3: best log2c:',num2str(bestLog2c),' best log2g:',num2str(bestLog2g),' accuracy/rmse:',num2str(bestcv),'%']);

% Plot the results
figure;
imagesc(cvMatrix); colormap('jet'); colorbar;
set(gca,'XTick',1:numLog2g)
set(gca,'XTickLabel',sprintf('%3.1f|',log2g_list))
xlabel('Log_2\gamma');
set(gca,'YTick',1:numLog2c)
set(gca,'YTickLabel',sprintf('%3.1f|',log2c_list))
ylabel('Log_2c');

end